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Abstract —We present a novel algorithm for reducing the 
state dimension, i.e. order, of linear parameter varying (LPV) 
discrete-time state-space (SS) models with affine dependence 
on the scheduling variable. The input-output behavior of the 
reduced order model approximates that of the original model. 
In fact, for input and scheduling sequences of a certain 
length, the input-output behaviors of the reduced and original 
model coincide. The proposed method can also be interpreted 
as a reachability and observability reduction (minimization) 
procedure for LPV-SS representations with affine dependence. 

I. INTRODUCTION 

In control applications, it is often desirable [16], [14] to 
use discrete-time linear parameter-varying state-space repre¬ 
sentations with affine dependence on parameters (abbreviated 
as LPV-SS representations in the sequel) of the form: 

x(t + 1)= A(p(t))x(t)+B{p(t))u(t) 
y{t)= C(p(t))x(t), 

where t G N, x{t ) G R ,!jc is the state, y(f) G R"? is the output, 
u(t) G R"“ is the input, and p(t) = \p\{t) ••• Pn p {t)] G 

P C R"p is the scheduling signal at time t G N. Here P is an 
arbitrary but fixed subset of IR"' 1 with a non-empty interior, 
and N denotes the set of natural numbers including zero. The 
matrices A{p{t)), B(p(t)), C{p{t)) in ([]} are assumed to be 
affine and static functions of p(t) of the form: 

n p 

A(p(t)) =A 0 + Y j A i p i (t), 

i= 1 
n p 

B(p(t))=B 0 + Y, B iPi(t), ( 2 ) 

t=l 

n p 

C{p(t))=Cv + Y d CiPi{t), 
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where A,- G R %x ”- C , B, G R' ljcX ”", C, G M. n y xn * are constant 
matrices for all i G {0,1,... ,n p }. 

Contribution of the paper Consider a LPV-SS represen¬ 
tation E of the form i|T]i and fix a positive integer N. In this 
paper, we present a procedure for computing another LPV-SS 
representation 

x(f + l) =A{p(t))x(t)+B{p(t))u(t) 
y{t)=C(p(t))x(t), 

such that for x(0) = 0, y(f) = y(t) for 0 < t < N, for all 
scheduling sequences (p(Q),p(l),...,p(N)) G P^ 1 and in¬ 
put sequences u = (m(0),m(1), ... ,u(N)) G (R' ! ") lv . Moreover, 
the state space dimension of £ is smaller than or equal to 
the state space dimension of E. In other words, given an 
LPV-SS representation E of order n x (state space dimension 
n x ) and a N G N\ {0}, we would like to find another LPV- 
SS representation E of order r < n x which has the same 
input-output behavior for all scheduling and input sequences 
of length up to N + lQ In addition, we would like the 
representation E to be a “good” approximation of E in 
terms of input-output behavior, even for scheduling and input 
sequences of length greater than N + 1 (see Remark □ lor 
what is meant by “good” here). Intuitively, it is clear that 
there is relationship between N and r. larger N yield a better 
approximation of the original input-output behavior, but they 
also result in larger values of r. In this paper, this relationship 
will be made more precise. Finally, by making use of this 
relation, the number N can be guaranteed to be chosen such 
that the resulting representation is a complete realization of 
the original model and it is reachable and/or observable. 
Therefore, the procedure stated in the present paper can also 
be used for reachability or observability reduction (hence, 
minimization) of an LPV-SS representation. 

Motivation LPV-SS representations are used in a wide va¬ 
riety of applications, see for instance [10], [19], [4], [18], [5]. 
Their popularity is due to their ability to capture nonlinear 
dynamics, while remaining simple enough to allow effective 
control synthesis, for example, by using optimal „ 

control. Model Predictive Control or PID approaches. LPV- 
SS representations arising in practice, especially which arise 
from first-principles based modeling methods, often have a 
large number of states. This is due to the inherent complexity 
of the physical process whose behavior the LPV-SS rep¬ 
resentations are supposed to capture. Unfortunately, due to 
memory limitations and numerical issues, the existing LPV 

1 Note that finding a representation £ with the same number of states as £ 
is in fact not necessarily useful, but it can happen that the proposed method 
does not allow us any other option. 


controller synthesis tools are not always capable of handling 
large state-space representations [8]. Moreover, even if the 
control synthesis is successful, large plant models lead to 
large controllers. In turn, large controllers are more difficult 
and costly to implement, and they often require application 
of reduction techniques. For this reason, model reduction of 
LPV-SS representations is extremely relevant for improving 
the applicability of LPV systems. 

To the best of our knowledge, the results of this paper 
are new. The tools which have been used in this paper 
stem from realization theory of LPV-SS representations [12], 
[17]. Similar tools were used for linear switched systems 
in [2]. In fact, we use the relationship between LPV-SS 
representations and linear switched systems derived in [12] 
to adapt the tools of [2] to LPV-SS representations. The 
method employed in this paper is related to that of [17]. The 
main difference is that [17] requires the explicit computation 
of Hankel matrices of LPV-SS representations. It should 
be noted that the size of the partial Hankel matrix of an 
LPV-SS representation increases exponentially (this will be 
stated more clearly in the paper, after necessary definitions 
are made). In contrast, the algorithm proposed in this paper 
does not require the explicit computation of Hankel matrices, 
and its worst-case computational complexity is polynomial. 
We present an example where the algorithm of [17] is not 
feasible due to the large size of the Hankel-matrix, while the 
algorithm of this paper works without problems. 

Regarding the literature, model reduction problem of LPV- 
SS representations was investigated in several papers [6], [7], 
[1], [21], [20], but except [20] they are only applicable to 
quadratic ally stable LPV systems. The method of [20] is 
applicable to quadratically stabilizable and detectable LPV- 
SS representations. In contrast, this paper does not impose 
any restrictions on the class of LPV-SS representations. In 
[15] joint reduction of the number of states and the number 
of scheduling parameters has been investigated. However, 
the method of [15] requires constructing the Hankel ma¬ 
trix explicitly. Hence, it suffers from the same curse of 
dimensionality as [17]. In addition, the system theoretic 
interpretation of the algorithm is less clear. To sum up, the 
main advantages of the proposed model reduction algorithm 
are the following: 

. it is applicable to arbitrary LPV-SS representations, 

• it has a clear system theoretic interpretation, 

• its computational (time and memory) complexity is 
polynomial in the number of states. 

The main disadvantage of the presented method is the lack of 
analytic error bounds. Note, however, that even for classical 
linear systems, there exists no analytical error bounds for 
model reduction algorithms which are based on moment 
matching. 

Outline: In Section [HI we present the formal definition 
and main properties of LPV-SS representations. In Section 
Hill we recall the concept of sub-Markov parameters for LPV- 
SS representations and give the precise problem statement. In 
Section [IV] we present the moment matching algorithm. In 


Section [V] the algorithm is illustrated on numerical examples 
and its performance is compared with the one of [17]. 

II. DISCRETE-TIME LPV-SS REPRESENTATIONS 

In this section, we present the formal definition of discrete¬ 
time LPV-SS representations and recall a number of relevant 
definitions. We follow the presentation of [12], 

In the sequel, we will use 

(n y ,n„,u v ,{(A,',B,',Ci)}"j! 0 ), (4) 

or simply E to denote a discrete-time LPV-SS representation 
of the form ([]}. In addition, we use to denote the set 
Ijj = {s G N | it < s < S 2 j-- An LPV-SS representation E is 
driven by the inputs {u(k)}°£ = 0 and the scheduling sequence 
In the sequel, regarding state trajectories, the 
initial state x(0) for an LPV-SS representation is taken to 
be zero unless stated otherwise. This assumption is made 
to simplify notation. Note that the results of the paper can 
easily be extended for the case of non-zero initial state. 

Notation 1: We will use II to denote the set of all maps 
of the form / : N —> H where H is a (possibly infinite) set. 
Using this, the sets //, fid, & and ffff are defined as °// = 
U N , 2d = P N , = T N and 32 = where U = IT", P = 
P C R»p, Y = JET-v and X = R n *. 

Consider an initial state xo € IR” J: of the LPV-SS repre¬ 
sentation E of the form (UJ. The input-to-state map AV t(1 : 
^ x 2d —> 32 and input-output map Yz. Xo '■ % x 2d —> <3d of 
E corresponding to this initial state xo are defined as follows: 
for all sequences u = {«(&)}>=() £ and p = {p(k)}^ =0 £ 
2*, let V£ rVo (u,p)(f) = x(f) and Fz lXo (u,p)(f) =y(t), t £ N, 
where x(t), y(t) satisfy ([]} and x(0) =xo . In the sequel, we 
will use Xv and IV to denote Xz.q and Yz.n respectively. That 
is, Xz and IV denote the input-to-state and input-output maps 
which are induced by the zero initial state. In fact, in the 
sequel we will be dealing with those input-output maps of 
LPV-SS representations which correspond to the zero initial 
state. 

The definition above implies that the potential input-output 
behavior of an LPV-SS representation can be formalized as 
a map 

(5) 

The value /(u,p)(?) represents the output of the underlying 
black-box system at time t, if the initial state x(0) = 0, 
the input u = {u(k)}^ =Q and the scheduling sequence p = 
{p(k)}~ =0 are fed to the system. Note that this black-box 
system may or may not admit a realization (description) by 
an LPV-SS representation, but the input-output behavior of 
any LPV-SS can be represented by a function of the form 
©. Next, we define when an LPV-SS representation realizes 
(describes) /. The LPV-SS representation E of the form (© is 
a realization of a map / of the form ©, if / equals the input- 
output map of E, i.e., / = IV. Two LPV-SS representations Ei 
and E 2 are said to be input-output equivalent if Yz l — Yz 2 . Let 
E be an LPV-SS representation of the form (|T}. We say that E 
is reachable, if R"' = span{AV(u,p)(f) | (u,p) £ x ffd,t £ 
N}, i.e. R' ijc is the smallest vector space containing all the 


states which are reachable fromx(O) = 0 by some scheduling 
sequence and input sequence at some time instance t, where 
ISN. We say that E is observable if for any two initial states 
x\,X 2 £ K"- v , — ^r,A' 2 implies jq = x 2 . That is, if any two 

distinct initial states of an observable E are chosen, then for 
some input and scheduling sequence, the resulting outputs 
will be different. 

Consider an LPV-SS representation Ei of the form (]T|i and 
an LPV-SS representation E 2 of the form 

E 2 = (n y ,n u ,n x ,{(A“,B“,C“)}Zo)- 

A nonsingular matrix y £ R"- rX ” A ' is said to be an LPV-SS 
isomorphism from Ei to E 2 , if for all i £ Iq P 

A“y = sfAi, b“ = y’Bi, c?y = c t . ( 6 ) 

In this case Ej and E 2 are called isomorphic LPV-SS 
representations. The order of E, denoted by dim(E) is the 
dimension of its state-space. That is, if E is of the form (|T|. 
then dim(E) = n x . Let / be an input-output map of the form 
©. An LPV-SS realization E is a minimal realization of /, 
if E is a realization of /, and for any LPV-SS representation 
£ which is also a realization of /, dim(E) < dim(E). We say 
that E is minimal , if E is a minimal realization of its own 
input-output map Yz- From [12], it follows that an LPV-SS 
representation E is minimal if and only if it is reachable and 
observable. In addition, if two minimal LPV-SS realizations 
are input-output equivalent, then they are isomorphic. Note 
that we defined minimality and input-output equivalence in 
terms of the input-output map induced by the zero initial 
state, hence we disregard autonomous dynamics. 

III. MODEL REDUCTION OF LPV-SS 
REPRESENTATIONS: PRELIMINARIES 

In this section, the sub-Markov parameters of a realizable 
input-output map / and its corresponding LPV-SS represen¬ 
tation E will be defined, and the moment matching problem 
for LPV-SS realizations will be stated formally. To this end, 
we recall the concepts of an infinite impulse response (HR) 
representation of an input-output map [17] and the concept 
of sub-Markov parameters. 

Consider an LPV-SS representation E of the form ©, 
and consider its input-output map / = IV. Recall from [17] 
that for any input sequence u = {u(k)}f (j and scheduling 
sequence p = {p(k)}f =0 , 

/(u,p)(f) = T £ (u,p)(f) = Y,( h ™op)(t)u(t-m) (7) 

m =0 

for all t £ N where 

(h 0 op)(t) = 0, (hiop)(t) = C(p(t))B(p(t - 1)), 

Vm > 1 : (h m op)(t) = (8) 

C(p(t))A(p(t - 1)) • • • A(p(t - m + 1 ))B(p(t-m)). 

The representation above is called the HR of / = Yz. From 
® and ©, it can be seen that the terms ( h m op)(t ), m > 0 


can be written as follows: 


(h 0 op)(t) = 0, 

rip rip 

(h\op)(t)= £ £ C q B qo p q (t)p qo (t - 1) 

q=0 <?o=0 

(h m op)(t) = 


E L 52 52 CqAj l ■ ■ 'Aj m _ l B qo p q j\...j m _ iqo 

q=0ji=0 j m -l=0qo=0 


( 9 ) 


where p 0 (k) = 1 for all k £ If, and p q j v ~j m _ m = Pq(t)p jl ( t - 
1 ) •'' Pj m -1 (t ~ m + 1 )Pq 0 (j ~m). 

Now we are ready to define the sub-Markov parameters 
of E. To this end, we introduce the symbol £ to denote the 
empty sequence of integers, i.e. £ will stand for a sequence of 
length zero and we denote by y(I 0 p ) the set {£}U{yi ■ • • j m \ 
m> l,ji,... ,j m £ Iq P } of all sequence of integers from Iq P , 
including the empty sequence. If s £ o5^(Iq P ), then |s| denotes 
the length of the sequence s. By convention, if .v = £, then 
|s| = 0. The coefficients 

r lq,q 0 (£)=CqBq 0 , 

Pq,qo O' 1 ' ■' jm) = CqAj 1 '' ' Aj m B qo , 

m> 1; q,j\,...,j m ,qo £ Iq P appearing in © are called the 
sub-Markov parameters of the LPV-SS representation E. In 
the sequel, the sub-Markov parameters T)^ 1/0 (s), q,qo £ Iq P , 
s £ oS^(Iq P ), |s| =m will be called sub-Markov parameters 
of E of length m. The intuition behind this terminology is as 
follows: the length of a sub-Markov parameter is determined 
by the number of Aj matrices which appear in (flOl) as factors. 

Note the sub-Markov parameters do not depend on the 
particular choice of an LPV-SS representation, but on the 
choice of the input-output map (provided that we fix an affine 
depency of the matrices of the LPV-SS representation on the 
scheduling variable). From [12] it follows that if E|, E 2 are 
two LPV-SS representations with static affine dependence 
on the scheduling variable, then their input-output maps are 
equal, if and only if their respective sub-Markov parameters 
are equal, i.e. Y Z{ = Yz 2 <V=> Vs £ y(f 0 r ) : V q , qo (s) = 
rjq, qo (s). Note also that another way to interpret the sub- 
Markov parameters is that they correspond to the derivatives 
of / with respect to the scheduling parameters. 

Example 1: Let E = (n y ,n u .n x , {(A;. B,,C,j}? =0 ) be an 
LPV-SS realization of the map / = Yz- Then the output of 
E due to the input u = {«(k)}“ =0 and scheduling sequence 
P = {p(k)}T=o at ti me f = 2 will be 

Yz( u,p)(2) =y( 2) = Y,(hiO P )(2)-u(2-i) 

i =0 

= 0 + (hi op) (2) • u (2 - 1) + (h 2 op) (2) • u(2 - 2) 

= C(p)B(p(t-\))u(\)+C(p)A(p(t-\))B(p(t-2))u(Q) 

2 2 

= E E C q B q 0 Pq( 2 )Pqo( l ) u ( l ) 

<7=0 <7o=0 
2 2 2 

+ EE £cAV«(2K(iK(0)«(0). 

q=0 71=0 90=0 


Recall that po(k) = 1 for all k £ Ik. In addition, observe 
from ©, that the output y(t), for t > 1 of an LPV-SS repre¬ 
sentation corresponding to an input sequence u = {«(&)}“_ 0 
and a scheduling sequence p = {p(k)}°[ {) is uniquely deter¬ 
mined by the sub-Markov parameters of length up to t — 1 
i.e., only the sub-Markov parameters of length up to t — 1 
appear in the output y(t) (see Example [[] for an illustration). 
Hence, if the sub-Markov parameters of length up to t — 1 
of two LPV-SS representations E and £ coincide, it means 
that E and £ will have the same input-output behavior up 
to time 1 for arbitrary input and scheduling sequences. This 
discussion is formalized below. 

Lemma 1 (I/O equivalence and sub-Markov parameters): 
For any LPV-SS representations E],E 2 , 

V(u,p) £ W x ^,k£ Iq : Y Zl (u,p)(k) =Y Zl (u,p)(k) 

if and only if 

Vs £ ^(Iq P ) , q, qo e Iq P , |s| < t - 1: r]fy 0 (s) = rj^ Q (s) 

This prompts us to introduce the following definition. 

Definition 1: Let E be an LPV-SS representation of the 
form 0. An LPV-SS representation £ of the form (0 is 
called a A'-partial realization of / = Yy, for some N £ N, if 

Vs£ J^(Io P ), 4 ,q , o £Iq P ,|s| <N: rj£ go (s) = Vq, 9o (s) 

That is, E is an N-partial realization of / = Te, if sub-Markov 
parameters of Yy and Yy up to length N are equal. In other 
words, £ is an A'-partial realization of Yy, if 

C q B % = C q B qo , \/q,q 0 £ Iq", 

CqAji ■ • 'Aj k Bq o = C q Aj l ■ ■ ■Aj k B qo , ,Vk £ Ij , 

fiq,q 0 Ji,...Jk&f 0 P - 

The problem of model reduction by moment matching for 
LPV-SS models can now be formulated as follows. 

Problem 1: Let E be an LPV-SS representation and let 
/ = Yz be its input-output map. Fix A' £ N. Find another 
LPV-SS realization £ such that dim(E) < dim(E) and £ is 
an A'-partial realization of / = Te- 

In order to explain the intuition behind this definition, we 
combine [13, Theorem 4] and [12] to derive the following. 

Corollary 1: Assume that E is a minimal realization of 
/ = Te and N is such that 2dim(E) — 1 < N. Then for any 
LPV-SS representation £ which is an A'-partial realization of 
/, £ is also a realization of / = Te and dim(E) < dim(£). 

Remark 1: Corollary Q] implies that there is a tradeoff 
between the choice of N and the order of E. Assume E is 
a minimal realization of / = Te- If N is chosen to be too 
high, namely if it is such that N > 2n x — 1, then it will not 
be possible to find an LPV-SS representation which is an 
A'-partial realization of / and whose order is lower than n x . 
In fact, if the model reduction procedure to be presented in 
the next section is used with any input N > 2 n x — 1, then 
the resulting LPV-SS representation £ will be a complete 
realization of / = Te. However, the order of £ will be the 
same as the order of E (provided that E is minimal). This 
relation between N and n x gives an a priori idea of how well 
the input-output map of £ approximates that of E. More 


specifically, we can expect the output error Te — Tg to be 
smaller when N is increased, as long as N < 2 n x — 1. This 
error will be zero for N > 2 n x — 1, since in this case £ will 
be a complete realization of Te. 


IV. MODEL REDUCTION OF LPV-SS 
REPRESENTATIONS 

In this section, first, the theorems which form the basis of 
the model reduction by moment matching will be presented. 
Then the algorithm itself will be stated. In the sequel, the 
image (column space) and kernel (null space) of a real matrix 
M is denoted by im(M) and ker(M) respectively. In addition, 
rank(M) is the dimension of im(M). We will start with 
presenting the following definitions for LPV-SS realizations 
of the form 0. 

Definition 2 (N-partial unobserx’ability space): The 
A'-partial unobservability space C'yfE) of E is defined 
inductively as follows: 

w) = n ker ( c ?)= 

<?eC 

^jv(E) = <?o(X) n n ker(^v_t(E)A ; ), N> 1. 


( 11 ) 


From [11], [12], it follows that E is observable if and only 
if <?jv(E) = {0} for all N >n x — 1. 

Definition 3 (N-partial reachability space): The A'- 
partial reachability space E) of E is defined inductively 
as follows: 


■^b(£) = span (J im(B ?0 ), 




^;v(E)=^o(£)+ E MAj&N-lP)), N>1. 


( 12 ) 




where the summation operator must be interpreted as the 
Minkowski sum. 

Again, from [11], [12] it follows that E is span-reachable if 
and only if dim(^)v(E)) = n x for all N > n x — 1. 

Remark 2: Let E be a LPV-SS representation of the form 
0. Recall from [17] the definition of the A'-step extended 
reachability matrix Ay and the definition of the A'-step 
extended observability matrix On of E. It is easy to see 
that ker(<9y) = &n( E) and im(A;v) = ^jv(E). Following 
[17] define Hankel matrix Hn.n of an LPV-SS representa¬ 
tion E as Hn,n = OnRn ■ Note that //y y is of dimension 
n y (n p + 1) { — + ’ ) xn„(w p +l)( — ) ■ i.e. it 

is exponential in N. Recall that [17] proposes a Kalman- 
Ho like realization algorithm based on the factorization 
of %jv for some N. The problem with this approach is 
that it involves explicit construction of Hankel matrices. 
Consequently, in the worst-case, memory-usage and time 
complexity of the algorithm [17] are exponential A'. In [17], 
N is chosen so that rank of //y y equals some integer n and 
the order of the LPV-SS computed from Hn.n will t> e at 
most n. While for many example, N will be small, it can 
happen that N is large, with N = n — 1 being the worst- 
case scenario, see Section [V] for an example. In addition. 




the method in [17] does not solve Problem [Q instead it 
relies on an approximation which is similar to balanced 
truncation. It yields an LPV-SS representation whose sub- 
Markov parameters are close to the corresponding sub- 
Markov parameters of the original LPV-SS representation. 
In Section [V] these remarks will be illustrated by numerical 
examples. 

Theorem 1: Let E = (n y ,n u ,n x , {(A ( -,B,-,C ; )} != ^ 0 ) be an 
LPV-SS representation, let V £ M n,;/r be a full column rank 
matrix such that 

E) = im(V). 

If E = (, n y ,n u ,r,{(Ai,Bj,Ci )}"= 0 ) is an LPV-SS representation 
such that for each i £ I 0 P , the matrices A ,C, are defined 
as 

A t = V~ l A t V, Bj = V~ l Bi, Q = C,V, 

where V 1 is a left inverse of V, then E is an /V-partial 
realization of the input-output map / = of E. 

This theorem follows from [2], [3] using [12]. For the sake 
of completeness, we present the proof below. 

Proof: Let V = 0. Since the conditions of Theorem Q] 
imply im(B ?0 ) C im(V), qo £ Iq P and V -1 is a left inverse 
of V, it is a routine exercise to see that VV~ l B qo = B qo . If 
V > 1, then im(A ;i . • • • Aj t B qo ) is also a subset of ^v(E) = 
im(V), i— 1 Hence, by induction we can show 

that VV~ x Aj j ---Aj l B clQ = Aj t - ■ ■Aj l B qo , i = 1, ...,1V, which 
ultimately yields 

V^j N A/, B qo = Aj N ■■■Aj l B qo . (13) 

Using (fl3l >. and C q = C q V, q £ Iq P , we conclude that for all 
i<N; q,q 0 ,ji,...,ji £ t Q p , 

Cqh • ■ ' Ajy B q Q = C q Aj i ■ ■ ■A j ] Bq (j 

from which the statement of the theorem follows. ■ 

Note that the number r is the number of columns in the 
full column rank matrix V, hence r < n x . This fact leads E to 
be of reduced order if V is sufficiently small, see Corollary 
[Q Using a dual argument, we can prove the following. 

Theorem 2: Let E = (n y ,n u ,n x ,{(Ai,Bj,Ci)} i f :0 ) be an 
LPV-SS representation, and let W £ W XHx be a full row rank 
matrix such that 

^v(E) =ker(W). 

Let W^ 1 be any right inverse of W and let 

T.= (n y ,n u ,r,{(A i ,Bi,Ci)} n i !L 0 ) 

be an LPV-SS representation such that for each i £ Iq P , the 
matrices are defined as 

Aj = WAiW~\ Bi = WBi, Q = QW~ l . 

Then E is an V-partial realization of the input-output map 

f = Yz of E. 

The proof is similar to that of Theorem D] 

Finally, by combining the proofs of Theorem Q] and 
Theorem [2] we can show the following. 

Theorem 3: Let E = (n y ,n u .n X) ((A,-. Tq) be an 

LPV-SS representation, and let V £ R"' xr and W £ IR' x ' ! - r 


be respectively full column rank and full row rank matrices 
such that 

&n( E) = im(V), ^v(E) = ker(W) and rank(WV) = r. 

If E = (n y ,.r, {(Aj.Bj.C ,)}" p ,) is an LPV-SS representation 
such that for each i £ Iq P , ,4,./),, C, are defined as 

A,- = WAjV(WV)~\ Bi = WBj , Q = CjV(W ) -1 

then E is a 2V-partial realization of the input-output map 
f = Y z of E. 

Note that having a 2V-partial realization as an approxi¬ 
mation realization would be more desirable than having an 
V-partial realization, since number of matched sub-Markov 
parameters would increase. However, it is only possible to 
get a 2V-partial realization for the original model E when 
the additional condition rank(V) = rank(W) = rank (ILL) = r 
is satisfied. 

Now, we will present an efficient algorithm of model 
reduction by moment matching, which computes either an 
V or 2V-partial realization E for an / which is realized by 
an LPV-SS representation E. First, we present algorithms for 
computing the subspaces E) and <C(v(E). To this end, 
we will use the following notation: if M is any real matrix, 
then denote by orth(M) the matrix U such that U is full 
column rank, im(t/) = im(M) and U T U =/. Note that U can 
easily be computed from M numerically, see for example the 
Matlab command orth. 

The algorithm for computing V £ R"' xr such that im(V) = 
,5?/v(E) is presented in Algorithm [[] below. 


Algorithm 1 Calculate a matrix representation of ^v(E), 
Inputs: ({Ai,Bi} ie1s np) and V 

Outputs: V £ R ,!jcXr such that rank(V) = r, im(V) =^v(E). 

V := U 0 , U 0 := orth [fi 0 ••• B„ p \. 

for k = 1... V do 

V := orth([V A 0 V A,V ••• A„ p V]) 

end for 
return V. 


By duality, we can use Algorithm |T] to compute a W £ 
W xn * such that ker(W) = &n{ E), see Algorithmic] 


Algorithm 2 Calculate a matrix representation of ^)v(E) 
Inputs: {Aj,Cj}. and V 

Output: W £ R rx ' !j , such that rank(VV') = r, and kerf IT) = 

W). _ 

Apply Algorithm D] with inputs ({Aj.Cj} ie f P ) to obtain 
a matrix V. 

return W = V T . 


Notice that the computational complexity of Algorithm [T] 
and Algorithm |2] is polynomial in V and n x , even though 
the spaces of E) (resp. (^fv(E)) are generated by im¬ 
ages (resp. kernels) of exponentially many matrices. Using 








Algorithm 3 Moment matching for LPV-SS representations 

Inputs: E = (ny,n I 4 ,w x ,{(A i ,S i ,C/)}"^ 0 ), Mode G {R,0,T} 
and N G N. 

Output: £= (ny,n u ,r,{(Ai,Bj,Ci)}y o ). 

Using Algorithm |T][2] compute matrices V and W such that 
V is full column rank, W is full row rank and im(V) = 
^v(E), ker(W) = ^v(E). 

if rank(V) = rank(W) = rank(WV) and Mode = T then 
Let r = ran ML) and 

Ai = WAjV{WV)~\ Q = QV(WV)~\ 

Bj = WBj. 

end if 

if Mode = R then 

Let r = rank(V), V 1 be a left inverse of V and set 
At = V~%V, Q = QV , Bi = V~ l Bj. 

end if 

if Mode = 0 then 

Let r = rank(W) and let W^ 1 be a right inverse of W. 
Set 

Ai = WA t W~\ Q = QW-\ Bi = WBi. 

end if 

return £ = {ny,n u ,r,{(A i ,B i ,Ci)} n j L 0 ). 


Algorithms [I] and [2] we can formulate a model reduction 
algorithm, see Algorithm 0 

Theorems |T| - [3] imply the correctness of Algorithm [3 

Corollary 2: Using the notation of Algorithm [3] the 
following holds: If railk(' V ) = ran MW) = rank(WV) and 
Mode = T, then Algorithm [3] returns a 2/V-partial realization 
of / = Yz (if Mode = T and the rank condition does not 
hold, the algorithm returns nothing). Otherwise, Algorithm 
[3] returns an /V-partial realization of / = IV. 

Note that even if the condition rank(V) = ran MW) = 
rank(WV) does not hold. Algorithm [3] can always be used 
for getting an /V-partial realization, by choosing Mode = 0 
or Mode = R. 

Remark 3 (Minimization of LPV-SS representations): 
From [12], it follows that if N > n x — 1 then 

oo 

^(E) = £ J’KE) = 

(=0 

span{A E (u,p)(f) | (u,p) G x S? f > 0}, 

^)v(E) = f] = 

i=0 

{x G R"* | Fx^(u,p)(t) = 0,V(u,p) Gf x ^,Vf > 0}. 

In other words, an LPV-SS representation E of the form ([!]) 
is reachable if and only if the dimension of its /V-partial 
reachability space E) is n x for all N >n x — 1, and E 
is observable if and only if the dimension of its /V-partial 
unobservability space ^)v(E) is 0 for all /V > n x I. In 
addition from [12], it follows that E is a minimal realization 


of its own input-output map if and only if E is reachable 
and observable. Hence, using this fact and [11], [17], it 
can be shown that Algorithm Q] can be used as an order 
minimization algorithm. That is. Algorithm [3] can be used 
consecutively with the inputs N >n x — 1, Mode = R (in 
this case, the resulting E will be reachable and it will be a 
realization of / = 1/) and N>n x —1, Mode = 0 (in this case, 
the resulting £ will be observable and it will be a realization 
of / = IV) for reachability and observability reduction for E, 
respectively. In turn, the resulting representation £ will be a 
minimal realization of / = }/. 

Remark 4 (Order r of the reduced representation): A 
disadvantage of the model reduction algorithm proposed by 
this paper is that the order of the reduced model produced by 
the method is unknown a priori. Namely, only the number 
N is chosen by the user as an input to the procedure, and 
the order of the reduced LPV-SS representation for this N 
is unknown beforehand. However, this issue can easily be 
solved by slightly modifying the method according to the 
concept of nice selections [3]. We omit the details of this 
approach due to lack of space. 

V. NUMERICAL EXAMPLES 

In this section, initially, the method stated in the present 
paper is applied to Example 4 in [17] and the result is 
compared with the one given in [17]. For this, both proce¬ 
dures are implemented in Matlab. The codes and the data 
used for both examples in this section are available from 
https://kom.aau.dk/~mertb/. 

First, the algorithm is applied to get a 3rd order approxi¬ 
mation to the LPV-SS realization of order 4 in Example 4, 
[17]. The original LPV-SS representation used in this case is 
of the form E = (n y ,n u ,n x , {(A,- ,B,-,C;)}”/ 0 ) with n y = n u = 1, 
n x = 4 and n p = 3. When N is chosen to be 1 and Mode = 
Reach, the resulting reduced order model £ is a 1-partial 
realization of 5/ of order 3. The scheduling signal used 
for simulation is of the form p{t) = [/3 y/—p sin(p)] 
where the parameter p takes its values randomly at each time 
instant, in the interval [—2^,0]. In addition, a white input 
u{t) ~ ,/U(0,1) is used. The upper limit of the simulation 
time interval is chosen to be N + 50 = 51. Since N = 1, 
the sub-Markov parameters of length at most 1 are matched 
with the original LPV-SS model E. The precise number of 
matched sub-Markov parameters is thus: 

(w P + 1 ) ( + ^ +1 ~ 1 ) (n P + 1 ) = 80 (14) 

The original model E and the the reduced order model £ 
are simulated for 500 different scheduling and input signal 
sequences of the type explained above, and their outputs y(t) 
and y(t) are compared for 1 =0,1,....//, where K is the 
number of steps of the simulation. For each simulation, the 
responses of E and £ are compared with the best fit rate 
(BFR) (see [9], [17]) which is defined as 

\l'Lf=o\\y( t )-ym\\ 2 2 j 


BFR = 100% max 









TABLE I 

Comparison of Alg.[3]and the Alg. in [17] 


The Proc. 

Mean BFR 

Best BFR 

Worst BFR 

Run Time 

Aiim 

76.5710% 

86.5821% 

64.9409% 

0.0430 s 

Alg. in [17] 

75.4364% 

85.4157% 

58.5798% 

0.0711 s 


Responses of the original and reduced order models 



Fig. 1. The responses of the original LPV-SS model E of order 4 and the 
reduced order LPV-SS model £ of order 3 acquired by Algorithm [3 The 
BFR for this simulation is = 76.5773%. 


where y m is the mean of {y(f)}£L 0 - 

For this example, the algorithms stated in this paper and 
in [17] are implemented for comparison. The mean of the 
BFRs, which is computed over 500 simulations, can be 
seen on Table U In addition the best and worst BFRs over 
500 simulations and the run-times for one single reduction 
algorithm are also shown in Table |I] The outputs y(t) and 
y(t) of the simulation which give the closest value to the 
mean of the BFRs are shown in Fig. Q] We used Algorithm 
[3]to perform model reduction using moment matching. From 
Table |T] it can be seen that both algorithms result in almost 
the same fit rates, whereas the algorithm stated in the present 
paper provides a 50% in terms of computational complexity. 

Next, a numerical example is presented to further illustrate 
the difference between the algorithms of the present paper, 
and the algorithm in [17]. The algorithm in the present 
paper is applied to get a reduced order approximation to a 
minimal LPV-SS model whose linear subsystems are stable. 
The original LPV-SS model used in this case is of the form 
E = {n y ,n u ,n x ,{(Ai,Bi,Ci)}y o ) with n y = n u = 1, n x = 7 and 
n„ = 5. The {(A ; -,B,-,C,)}T 0 parameters of E are as follows: 


A 0 = 

-0.5 

0 6 xl 

0.5471 

06x1 

0lx5 

06x5 



Olxl 

Olxl 

Olxl 

0lx4 

Ai = 

Olxl 

0.3 

0.2285 

0lx4 


0 5 xl 

O 

L/t 

X 

O 

Lft 

X 

05x4 


0 2 x 1 

X 

<N 

o 

0 2 xl 

0 2x 3 

a 2 = 

0lx2 

-0.4 

0.4741 

0lx3 


04x2 

0-lxl 

04x1 

04x3 


03x3 

o 

u> 

X 

03x1 

03x2 

0lx3 

-0.7 

0.9362 

0lx2 

°3x3 

03x1 

03x1 

03x2 

04x4 

o 

-L* 

X 

04x1 

X 

■'d- 

O 

0lx4 

0.5 

0.4367 

Olxl 

0 2x 4 

0 2x i 

0 2x l 

0 2 xl 

05x5 

X 

m 

O 

l 

X 

n 

O 


0lx5 

0.1 

0.0573 


0lx5 

Olxl 

Olxl 



Bo = [0 0 0 

B\ = [0 0 0 

b 2 =[o 0 0 

B 3 = [1 0 0 

B 4 = [0 1 0 

B S =[ 0 0 1 

Ci =[100 
where 0 aX b, a,b £ N\{0} 


0 0 0 1] T 
0 0 1 o] T 
0 1 0 o] T 
0 0 o o] T 
0 0 o o] T 
0 0 o o] T 

0 0 0 0] ,Vi £ Iq p . 
denotes the zero matrix of dimen¬ 


sion ax b. 

The resulting reduced order model E is a 2-partial realiza¬ 
tion (hence N = 2) of TV of order 3. A random scheduling 
signal and u(t) ~ (0,1) is used for simulation. The upper 

limit of the simulation time interval is chosen to be N + 50 = 
52. Since N = 2, the sub-Markov parameters of length at 
most 2 are matched with the original LPV-SS model E. Note 
that the precise number of matched sub-Markov parameters 
can be found by using (fT4t with n p = 5, N = 2, which is in 
this case 1548. 

The output y(t) of the original model E and the output 
y(t) of the reduced order model £ are simulated again for 
500 random scheduling and white Gaussian input signal 
sequences. For this example, the mean of the BFRs over 500 
simulations is 93.4888%; whereas, the best BFR is 99.3192% 
and the worst is 47.9013%. The elapsed time for all of the 
simulations is 5.347983 seconds. The outputs y(t) and y(t) 
of the simulation which gives the closest value to the mean 
of the BFRs are shown in Fig. [2] It can be seen that the 
responses of both models are exactly matched until (and 
including) the time instant 1 = N = 2. In addition. Fig. [2] 
together with its BFR= 93.4993% show that the reduced 
order model E captures the behavior of the original model 
E accurately, for the rest of the total time horizon, i.e., for 
t > 2. Note that when the method given in [17] is applied 
to this example, the method breaks down by running out of 
memory while trying to compute the smallest Hankel matrix 
of rank dim(E) =n x =l. 

For the same example, the procedure in the present paper 
is applied again with N = 4. Then the resulting reduced order 
model is of order 5 and the BFRs over 500 simulations 
are as follows: Mean BFR= 97.4010, Best BFR= 99.8494, 
Worst BFR= 75.9829. The elapsed time for all simulations is 
5.621238 seconds. The outputs y(t) and y(t) of the original 
and the reduced order models which give the closest value 
to the mean of the BFRs are shown in Fig. 0 Finally, 
























Responses of the original and reduced order models 



Fig. 2. The responses of the original LPV-SS model E of order 7 and the 
reduced order LPV-SS model E of order 3 acquired by Algorithm [3] The 
BFR for this simulation is =93.4993%. 


Responses of the original and reduced order models 



Fig. 3. The responses of the original LPV-SS model £ of order 7 and the 
reduced order LPV-SS model £ of order 5 acquired by Algorithm [3] The 
BFR for this simulation is = 97.4002%. 

the procedure is applied to get a full realization. For this 
example, for N >6, the reduced LPV-SS representation has 
the same order with and it is isomorphic to the original LPV- 
SS representation considered. Hence, it is a full realization 
of / = Yv. The elapsed time for computing one such full 
realization for this example is 0.029692 seconds Note that 
this is the run-time for only one reduction procedure. No 
simulations were done to compare the outputs in this case, 
because they would be exactly the same for all input and 
scheduling sequences. 

VI. CONCLUSIONS 

A model reduction method is presented for discrete time 
LPV-SS representations with affine static dependence on the 
scheduling variable. The method makes it possible to find a 
reduced order approximation to the original LPV-SS model, 
which has the same input-output behavior for scheduling 
and input sequences of a pre-defined, limited length. The 
presented method can also be used for reachability and ob¬ 
servability reduction (i.e., minimization) for LPV-SS models. 
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